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q , 1 Core-collapse Supernovae: A Numerical Challenge 

(N : 

The approach to radiative transfer described in this contribution is being 
y—i • developed for use in simulations of core-collapse supernovae. These events are 

the deaths of stars more than about eight times as massive as the Sun, caused 
by the catastrophic collapse of the star's core. This collapse is triggered in 
part by electron capture on heavy nuclei, resulting in the emission of weakly- 
■^j- ■ interacting particles called "neutrinos." Early in the collapse process these 

| neutrinos escape freely, but eventually densities are sufficiently high that even 

neutrinos are trapped. Collapse is finally halted when central densities reach a 
few times the matter density of atomic nuclei. The newly-born neutron star — 
• the compact object resulting from this process — is a hot thermal bath of dense 

nuclear matter, electrons, positrons, neutrinos, and antineutrinos. Neutrinos 
Q ■ and antineutrinos, having the weakest interactions among the species present, 

^—j | are the most efficient means of cooling; their emission accounts for virtually 

all of the gravitational potential energy released during collapse. 
| Because neutrinos dominate the energetics of the supernova process, neu- 

trino radiative transfer is a central feature of the core collapse phenomenon. In 
k>( | particular, energy transfer from neutrino radiation to infalling stellar matter 

■ may be crucial to the supernova explosion mechanism ^ [2] . In addition to 

energy, neutrinos exchange a quantity called "lepton number" with the fluid; 



this composition variable affects the fluid's equation of state, and also has a 
strong influence on the relative abundances of nuclear species synthesized in 
the supernova environment. 

The radiative transfer of energy and lepton number is of particular interest 
in the semi-transparent region near the neutron star surface. This is a region 
of transition from the optically thick interior — where the diffusive neutrino 
field is nearly isotropic — to the optically thin exterior, where the neutrino ra- 
diation becomes strongly forward-peaked. Hence energy- and angle-dependent 
neutrino transport is key to accurate modeling of core-collapse supernovae. 
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At least three groups have published reports of spherically symmetric core- 
collapse supernova simulations with energy- and angle-dependent neutrino 
transport. Two of these groups — centered at the Max Planck Institute for 
Astrophysics |3] and the University of Arizona 0] — employed a method in 
which angle-integrated and angle-dependent neutrino transport equations are 
iterated to simulataneous convergence. 3 In contrast, the group centered at 
Oak Ridge National Laboratory 5 performed a direct solution of the angle- 
dependent transport equation, fully discretized in all variables in all terms. All 
these efforts were grid-based (with the first two groups employing "tangent 
rays" to discretize the momentum angle). 

The iteration of angle-integrated and angle-dependent transport equations — 
which might be called an "iterated moment method" — works as follows. Ze- 
roth and first angular moments of the neutrino transport equations are formed, 
with energy dependence retained. The (energy-dependent) zeroth and first 
angular moments of the neutrino distribution thus become variables to be 
evolved. The second and third angular moments also appear in these equa- 
tions; to close the system at the first moment, the higher angular moments 
are expressed as numerical factors (the so-called "Eddington factors") mul- 
tiplying the zeroth moment. The Eddington factors can be computed from 
the solution of the angle-dependent transport equation; this is solved with 
a simplified collision integral, which is expressed in terms of the zeroth and 
first angular moments of the neutrino distribution. In summary: The mo- 
ment equations need Eddington factors for closure, which are obtained from 
the solution of the (simplified) angle-dependent transport equation; while the 
angle-dependent transport equation requires the zeroth and first moments for 
its simplified collision integral. This system is iterated to convergence. 

The "direct method" of solving the transport equation is not subject to a 
structural limitation of the iterated moment method. In the direct method, 
all terms are discretized in all variables — time, space, energy, and angles. In 
particular the angle dependencies of neutrino scattering and pair production 
terms are fully represented, while in the iterated moment method only the 
I = 0, 1 terms in a Legendre expansion of these collision kernels are em- 
ployed; these are the only terms in the expansion that can be constructed 
from the zeroth and first moments of the neutrino distribution. (In principle, 
the Eddington factors for the second and third moments could be used to 
get two additional terms in the Legendre expansions, but this has not been 
implemented to date.) In the supernova environment, these truncated angular 
expansions might not be accurate representations of neutrino pair production 
and neutrino scattering by electrons and positrons — important processes in 
determining the spectra of neutrinos emerging from the nascent neutron star. 
(Another approximation of undemonstrated safety — simplification or neglect 
of angular aberration in the angle-dependent transport equation [31 El — seems 
less fundamental to the iterated moment method itself.) 



3 See also H.T. Janka's contribution to this volume. 



An Approach to Neutrino Radiative Transfer in Supernova Simulations 



3 



While energy- and angle-dependent neutrino transport is an important 
advance even in spherical symmetry, the physics of stellar core collapse de- 
mands a spatially multidimensional treatment. There is ample theoretical and 
observational evidence for this conclusion (see for example 6 ). 

Inclusion of energy- and angle-dependent neutrino transport in spatially 
multidimensional simulations represents a significant computational challenge. 
Consider the "direct method" mentioned above, in which the neutrino distri- 
bution functions and all terms in the transport equation are discretized in 
all variables. Assume azimuthal symmetry; let the numbers of spatial zones 
in (r, 6) be (256,128), and the numbers of momentum bins in energy and 
angle variables (e,#, (p) be (64,32,16). An implicit time evolution algorithm 
requires the solution of a large linear system, in which the matrix represents 
the coupling between values of the neutrino distribution function at different 
points on this five-dimensional grid. The spatial coupling of points having the 
same momenta is very sparse, involving only nearest neighbors. But neutrino 
scattering and pair production terms in the collision integral involve dense 
and extended coupling in momentum space; this means that the matrix con- 
tains 256 x 128 dense blocks of size (64 x 32 x 16) 2 . Storing all of these dense 
blocks in double precision would require a few 10 14 bytes, beyond the capacity 
of today's terascale machines. Moreover, the elements of the dense blocks are 
expensive to compute, rendering unattractive an approach in which they are 
generated "on the fly" several times during each solution of the linear system. 
These considerations motivate the algorithm presented in this paper. 

2 Formalism for Radiative Transfer 

Some details of a radiative transfer formalism and its finite differenced rep- 
resentation will now be presented. Conservative formulations of radiative 
transfer — motivated by the importance of accurately tracking energy and lep- 
ton number transfer in the supernova environment — are discussed first. This 
is followed by a discussion of the finite differencing of the terms that do not 
depend on the velocity of the background fluid. 

2.1 Conservative Formulations of Relativistic Radiative Transfer 

A radiative transfer calculation involves a variable describing the distribution 
in phase space of the particles comprising the radiation. 

The particle distributions we consider in some detail are the scalar distri- 
bution function / and the specific number density J\f; we also mention the 
specific energy density T. Each of these is taken to be a function of spacetime 
coordinates and three- momentum variables u l . (Greek and latin letters 
are spacetime and space indices respectively. Hatted indices indicate quanti- 
ties measured in an orthonormal frame comoving with the fluid with which 
the particle species interact, and unadorned indices indicate components with 
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respect to a global coordinate basis.) The momentum variables u 1 arise from 
a change of variables (e.g. to momentum space spherical coordinates) from 
p l , the Cartesian spatial momentum components measured in a comoving 
orthonormal frame. 

The scalar distribution function / gives the the number of particles dN in 
an invariant spacetime 3-volume element dV and invariant momentum space 
volume element dP jZj: 

dN = f(x,p) {-v -p)dVdP. (1) 

The quantities x, v, and p are 4- vectors, and p is the spatial 3-vector portion 
of p. The unit 4- vector v is timelike, and defines the orientation of dV : 

dV — s/^ge^vpaV 11 d\x v d2X p d^x 17 , (2) 

where g is the determinant of the metric tensor (taken to have signature 
— h ++) and e.p Vpa is the Levi-Civita alternating symbol (eoi23 = +1). The 
momentum space volume element is 

, D , d 1 p l d 2 p j d 3 p k 

(-Po) 

d^du^du 3 , (3) 



det 



dp 
du 



where specialization to momentum variables u l in the comoving frame has 
been made in the second line, and 



e = p° = V|P| 2 +™ 2 (4) 

is the particle energy measured in the comoving frame (m is the particle mass). 
The definition of / in eq. makes a convenient connection to nonrelativistic 
definitions of the distribution function; for an equivalent but more geometric 
approach see Ref. |S] ■ The transport equation for electrically neutral particles 
satisfied by / is the Boltzmann equation [71 |S1 OH E| : 

In this expression, is the transformation between the coordinate basis 
and the comoving orthonormal basis; it involves a transformation to an or- 
thonormal "lab frame" basis followed by a Lorentz boost to the orthonormal 
comoving frame. The connection coefficients in the orthonormal comoving 
basis are 

± vp — fj,*^ v p J up \ p*^ p Qj,p ' ^ ' 

Because the transformation D 1 ^ contains a Lorentz boost, the term dD 1 p/dx p 
gives rise to Doppler shifts and angular aberrations associated with this trans- 
formation. The coordinate basis connection coefficients, 
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1 vp ~ 2 9 



dg vp 

V dxP ' dx v dx a 



(7) 



where g^ are the metric components, give rise to energy shifts and angular 
aberrations associated with spacetime curvature and the use of curvilinear 
coordinates. The right-hand side of JHJ is the invariant collision integral. 

The Boltzmann equation can be put in a number-conservative form 10 , 
which motivates the definition of the specific particle number density TV. 
Specifically, eq. (JSJ can be rewritten as 



d 



-g dx^ 



<l,i|£ 

du 




t 5 -r>v> 



<1...|£ 

du 




= C[/]. 



(8) 



This form is called "conservative" because the left hand side is expressed as 
divergences in spacetime and momentum space, so that volume integrals of 
these terms are transparently related to surface terms. In particular, it is 
obvious that the momentum space divergence vanishes upon integration over 
dP (given by eq. ©), yielding the equation for particle number balance: 



1 



d 



-g dx^ 



(V=?jv) = / c[/]dP, 



where 



ATM 



(9) 



(10) 



are the coordinate basis components of the particle number flux vector. This 
motivates the definition of the specific particle number density Af, given by 



Af: 



P 



AC , 



(11) 



From eqs. (|10|) and l|ll|). we see that the specific number density is the con- 
tribution of each comoving frame momentum bin to the lab frame number 
density: 

iV° = [ AT edP. (12) 



In terms of Af, eq. JHJ can be rewritten as 

1 d 



<l,l [ £ 
du 




-g dx^ 

det[^ 
du 



^H c - ,i3) 



where we have defined 
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Note that is not a 4-vector. 

Another reformulation of the Boltzmann equation is an energy-conservative 
form JU| , which motivates the definition of the specific particle energy density 
T. This reformulation is not detailed here; suffice it to say that its momen- 
tum integral gives a transparent connection to the divergence of the particle 
stress-energy tensor. 

The number- and energy-conservative formulations — which respectively 
facilitate computation of lepton number and energy transfer to essentially 
machine accuracy — might be used in a couple of different ways. Both formu- 
lations could be solved separately, in order to nail down the transfer of both 
energy and lepton number. The values of the scalar distribution function im- 
plied by these two different distributions would then serve as a consistency 
check. Alternatively, only one of the formulations might be solved, with the 
analytic relationship between the two conservative formulations 10 being 
used to design finite difference expressions that provide consistency with the 
other formulation. This latter philosophy has been employed in the work of 
the group centered at Oak Ridge 5 . 



2.2 Finite differencing of selected terms 

Consider I jl^jl in two spatial dimensions in spherical coordinates, with the 
assumption of a static (zero velocity) background and flat spacetime: 

dJV cos $ d , 9 , rt sin d cos ip d , . n . _ 

-qT + -^Tjr (r AO + ^7^™ (sin AO - 

ot r z or ' r sin oO 

—L-JL (sin 2 -— ^-(sintfsin^AT) = I C . (15) 

In this equation, (r, 9) are the spatial radius and polar angle, ip) are mo- 
mentum space angles, and e is the particle energy. 

There are two things to keep in mind in constructing a finite-differenced 
representation of (|f 5(1 . First, one can take advantage of the conservative form 
to make numerical "volume integrals" transparently related to numerical "sur- 
face integrals." Second, notice that for M spatially homogeneous, the second 
and fourth terms of l|I5l) cancel, as do the third and fifth terms. The finite 
difference representation should respect this cancellation. 

A conservative differencing of the spatial divergence in (|I5|I is 

L -(cos 0V [{A r M) i+ x,j> - (A r J\f) itj >] + 



Vi 



' (sin #)/3' (cos y>) 7 ' [(Ae N)i> - {AeN)i>,j\ , (16) 



where the geometric factors are 
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V V j, =2ir(r l ,) 2 sm{6 r )(Ar) l ,(A0) J ,, (17) 
(A-)ij' = Mnf sin(^)(z}%, (18) 
{A e ) v ,j = 27rry smi^iAr)^. (19) 

In the unite-difference expressions in this subsection, the subscripts index 
spatial zones in the (r, 9) spatial coordinate directions, and subscripts (a, f3, 7) 
index momentum bins in the (e, momentum space coordinate directions. 
Unprimed indices denote values evaluated on the surfaces ( "edges" ) of spatial 
zones or momentum bins. The values of the coordinates (r, 6, ip) on the 
zone and bin edges are prescribed by the user. Primed indices denote values 
evaluated at zone or bin centers. Given the "edge values," the coordinates of 
the zone centers are taken to be 



1 (Ji , „3 \V3 



— \ Ta +r 



2 

= arccos 



1+1) ' 
-(cos(0,-) +cos(0 i+1 )) 



(20) 
(21) 



and the zone widths are (Ar)? = rj+i — 7^ and (A9)jt = 9j+i — 0j. The 
definitions of (cost?)/?', (sin??)^/, and (cosy) 7 ' in 116(1 will be given below. 
Finally, the values of j\f on zone surfaces in (|16f) are given by a particular linear 
interpolation of zone center values. This linear interpolation — which depends 
on the neutrino mean free paths — has the effect of shifting from "diamond" 
differencing in diffusive regimes to "upwind" (or "donor-cell" ) differencing in 
free-streaming regions; see jS] for details. 

A conservative differencing of the momentum space divergence in (|15fl is 



1 f 1 



1 (cot 9)? 



(A$ AO/3+1,7' - (AjAO/5,7'] ~ 

{A v N)f3> n +i-{A v N) , (22) 



Va',0',7' T V 

where the momentum space "geometric" factors are 

V Q ^, 7 < = (e a ,) 2 S in(dp,)(Ae) a ,(A-&)p,(A<p) Y , (23) 
= (e a/ ) 2 [sm^ p )} 2 (Ae) a ,(A l p) Y , (24) 
(Ap)p' t i = (e a/ ) 2 sm(dp,)sm(ip 1 )(Ae) al (A^) l3l . (25) 

The bin center values e a i and dpi are given just as ry and 9j> in 1)20(1 and l|21|l. 
Also, the momentum bin widths are given by the difference of the bounding 
edge values, just as in the case of the spatial zones. The definitions of (l/r)i' 
and (cot 0)ji will be given below. The values of N on momentum bin edges 
are given by "upwind" or "donor-cell" differencing (see 

Finally, we show the remaining definitions needed to ensure that the sec- 
ond and fourth terms in (|15|l . as well as the third and fifth terms, cancel 
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for TV spatially homogeneous. Given the differencings defined above, this is 
accomplished with the following definitions: 

(sin%, =sin(iV), (26) 
/ <n sin(^ /3+1 ) 2 - sin(^) 2 

(C0S ^'= 2«nVr)(A*)r ' (27) 



1\ (n + i) 2 - (n) 2 



(29) 



3 Radiative transfer algorithm and distributed-memory 
implementation 

In this section classes of operators involved in radiation transport arc identi- 
fied, and a strategy for implementing them on distributed-memory computer 
architectures is described. 

The equations of neutrino radiative transfer are the transport equation 
for whatever radiation particle distribution function is used, together with 
equations that describe lepton number and energy transfer (the latter are 
given by appropriate momentum integrals of the transport equation). The 
terms in these equations correspond to operators acting on the discretized 
distribution function and transfer quantities. This can be expressed as 

F[y] = o, (31) 

where y denotes the set of unknowns at a given time step: .ZVspecies x -/V spacc x 
^momentum unknown values of the distribution functions of iV spe cies neutrino 
species in N spacc spatial zones and iV momcn t um momentum bins, and 2 x N space 
unknown values of energy and lepton number transferred to the fluid in the 
N space spatial zones. The total operator F has various pieces: 

F = T + S + M + C. (32) 

The time derivative operator T relates unknowns at fixed position x and mo- 
mentum u at different times t. The space derivative operator S is linear, and 
connects nearest neighbors in x at fixed u and t; similarly, the momentum 
derivative operator M is linear and connects nearest neighbors in u at fixed x 
and t. (The operators S and M are divergences in the conservative formula- 
tions.) In the case of astrophysical neutrino transport, the collision operator 
C is nonlinear due to neutrino-neutrino interactions and phase space block- 
ing associated with the Pauli exclusion principle; and because of scattering 
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and pair production and annihilation processes, it exhibits extensive, nonlocal 
coupling in u at fixed x and t. 

The large disparity between hydrodynamic (~ s) and neutrino inter- 
action (~ 10 -10 s) time scales in the collapsed core of the supernova environ- 
ment calls for implicit evolution of the transfer equations. This means that in 
a time step in which the system is evolved from time t n to t n+ , the operators 
S, M, and C are evaluated at t n+ . The discretized transfer equations are 
then 

T{y n+ \ y n ) + S(y n+1 ) + M{y n+1 ) + C{y n+1 ) = 0, (33) 

where we have used the notation y n = y(t n ) and suppressed the dependence 
on space and momentum variables. The dependence of T on the values of y 
at only two different times indicates that a method that is first order in time 
is being used (specifically, backward Euler). 

Because of the nonlinearity of the collision operator C, (|33|l is a set of 
nonlinear algebraic equations for the discretized values of the distribution 
functions and transfer variables; this system is solved with a Newton-Raphson 
iteration procedure. Specifically, l|3f (I is linearized: 

J ■ Ay = -F, (34) 

where J = dF/dy is the Jacobian matrix. In this linearized equation, J and 
— F are evaluated at a guess (y n+1 ) guess for the value of the distribution 
function at the new time t n+1 . The solution Ay of this linear system provides 
a new guess (y n+1 ) nC w guess = (j/" +1 )gucss + Ay. This procedure is iterated to 
convergence of y n+1 to the solution of |3*3jl. 

To solve the linear problem at the heart of each Newton-Raphson iteration, 
a simple fixed-point method is employed. The basic idea is as follows. Start 
with a guess (Z\y)g Ue ss for the solution of l|34|) . and compute the residual r, 

r=(-F)-J-(Ay) gucss . (35) 

Given (</ _1 ) a pprox, an approximate inverse of J, compute the correction c, 

C = {'J )approx ' ^- (36) 

Why the name "correction" ? Note that if one computes c oxac t with the exact 
inverse J -1 , lp?H|. and (JB3J give 

^cxact J ' ( F^ J J (Z\7/)g Uess 

= Ay - (Ay) sueas , (37) 

the difference between the exact solution Ay to eq. Ij34(l and (Ay) gucss . When 
c is computed with an approximate inverse it does not provide the exact 
difference between Ay and (Ay) guess , but it does provide a new (and hopefully 
improved) guess (Ay) ncw guoss = (Ay) gucss + c. This procedure is iterated to 
convergence of Ay to the solution of (|34(l . 
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Details of the implementation of this iterative fixed-point method for the 
solution of the large linear system will now be given. In practice, the inverses 
of two different approximate Jacobian matrices are applied in succession in 
each iteration (this also requires the computation of two residuals in each 
iteration). The first approximate Jacobian is 

'/momentum = Jt + Jm + JCi (38) 

which consists of the contributions to the Jacobian from the operators T, M, 
and C in l|33fl . As previously described, this combination of operators densely 
couples different momenta u at fixed spatial position x; hence ./momentum 
consists of iVspaco independent dense blocks — one for each spatial zone — and 
(^momentum) 1 consists of individual inverses of these dense blocks. With the 
spatial grid partitioned among the many processors of a distributed-memory 
computer, the inversion of these separate blocks is trivially parallelized. The 
second approximate Jacobian is 

■/space = Jt + Jst (39) 

arising from the operators T, S in (|33l) . By reasoning similar to above, J spa ce 
can be conceptualized as iV momontum independent matrices, but this time with 
sparse coupling, because the derivatives in S only require nearest neighbors 
in space. Having chosen to partition the spatial grid, parallel solution of these 
independent "spatial matrices" requires an "all-to-all" communication to give 
each processor all the spatial data for its share of momentum bins; but the 
fact the matrices are sparse makes this communication manageable. 

In addition to its simple structure, this fixed-point method has an impor- 
tant practical advantage over other iterative linear solver algorithms. As men- 
tioned at the end of the first section, in spatially multidimensional problems 
simultaneous storage of all the dense blocks is impractical. The fixed-point 
algorithm outlined above can be structured so that each processor can con- 
struct a few dense blocks at a time, use them in all steps required in a given 
iteration, and discard them. In contrast, other linear solver algorithms seem 
to require dense blocks to be discarded and rebuilt multiple times in each 
iteration. 

A code that implements the algorithm described above — written in For- 
tran 90, and using the MPI library for message passing — is being developed 
and tested at Oak Ridge National Laboratory, for eventual use in core-collapse 
supernova simulations. The implementation has been tested in both one and 
two spatial dimensions on a "homogeneous sphere" problem which has an 
analytic solution, with good results. The test has also been generalized to an 
"inhomogeneous sphere" problem, with emissivity and opacity varying in spa- 
tial polar angle. With regard to performance, it is found that inversion of the 
dense blocks dominates the computation; communication costs are not exces- 
sive. Because dense matrix solvers (e.g. the LAPACK library) are typically 
highly optimized, the dominance of the computation by dense blocks ensures 
that computational resources are used efficiently. 
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